#include "springSmooth.h"
#include "mesh.h"
#include "distance.h"
#include <math.h>
#include <fstream>
#include <iostream>
#include <iomanip>

springSmooth::springSmooth() : _k_frac(0.5), _nt(1), _t(0) {}

springSmooth::springSmooth(const bool &enableTimeStep, const double &residual_rel, const int &nt, const double &k_frac, const double &mag, const double &t) : _enableTimeStep(enableTimeStep), _residual_max(residual_rel), _nt(nt), _k_frac(k_frac), _mag(mag), _t(t)
{
}

springSmooth::~springSmooth() {}

void springSmooth::generateMesh()
{
    remove("tec.plt");
    grid.generateMesh();
}

void springSmooth::saveMesh()
{
    grid.saveMesh();
}

void springSmooth::getTopology()
{
    grid.getNeighborElem();
    grid.getNeighborNodes();
}

void springSmooth::run()
{
    grid.getTop();
    grid.getBottom();
    grid.getLeft();
    grid.getRight();

    // 获取数据
    unsigned int _points_num = grid.returnPointsNum();
    std::vector<std::vector<int>> neighborNode = grid.returnNeighborNodes();
    std::vector<std::vector<double>> mesh_points = grid.retrunMesh_points();

    std::vector<unsigned int> bottom_boundary = grid.return_bottom_boundary();
    std::vector<unsigned int> top_boundary = grid.return_top_boundary();
    std::vector<unsigned int> left_boundary = grid.return_left_boundary();
    std::vector<unsigned int> right_boundary = grid.return_right_boundary();

    d_previous.resize(_points_num);
    d_current.resize(_points_num);
    for (int i = 0; i < _points_num; i++)
    {
        d_previous[i].resize(2);
        d_current[i].resize(2);
    }
    _dt = _t / _nt;
    this->outNewMesh(0.0);
    // 设置边界位移
    for (int i = 0; i < bottom_boundary.size(); i++)
    {
        d_previous[bottom_boundary[i]][0] = 0;
        d_current[bottom_boundary[i]][0] = 0;
        d_previous[bottom_boundary[i]][1] = (_enableTimeStep ? _mag / _nt : _mag) * sin(mesh_points[bottom_boundary[i]][0]);
        d_current[bottom_boundary[i]][1] = (_enableTimeStep ? _mag / _nt : _mag) * sin(mesh_points[bottom_boundary[i]][0]);
    }

    std::vector<double> numerator = {0.0, 0.0};
    double kij = 0, denominator = 0;
    int counter = 0;
    double current_t = 0.0;
    for (int c = 0; c < _nt; c++)
    {
        current_t = (c + 1) * _dt;
        while (true)
        {
            for (unsigned int i = 0; i < _points_num; i++)
            {
                if (myfind(bottom_boundary, i) || myfind(top_boundary, i) || myfind(left_boundary, i) || myfind(right_boundary, i)) // 如果是bottom
                    continue;

                for (int j = 0; j < 10; j++)
                {
                    if (neighborNode[i][j] == -1) // 10 5 7 -1 -1 -1
                        break;
                    kij = _k_frac / sqrt(distance(mesh_points[i], mesh_points[neighborNode[i][j]])); // 这块浪费了计算效率
                    // 计算分子
                    for (int k = 0; k < 2; k++)
                    {
                        numerator[k] += kij * d_previous[neighborNode[i][j]][k];
                    }
                    denominator += kij;
                    // 计算分母
                }
                for (int k = 0; k < 2; k++)
                {
                    d_current[i][k] = numerator[k] / denominator;
                }
                std::fill(numerator.begin(), numerator.end(), 0.0);
                denominator = 0.0;
            }
            residual = calResidual(d_previous, d_current);
            d_previous = d_current;
            std::cout << "step" << c << ":" << residual << std::endl;
            if (counter == max_run_count || residual <= _residual_max) //
                break;
            counter++;
        }
        std::cout << "time step:" << c << "finished!" << std::endl;
        this->outNewMesh(current_t);
        // 从原始位置重新求解新的位移
        counter = 0;
        numerator = {0.0, 0.0};
        for (int i = 0; i < bottom_boundary.size(); i++)
        {
            d_previous[bottom_boundary[i]][0] = 0;
            d_current[bottom_boundary[i]][0] = 0;
            d_previous[bottom_boundary[i]][1] = (_enableTimeStep ? (c + 2) * _mag / _nt : _mag) * sin(mesh_points[bottom_boundary[i]][0]);
            d_current[bottom_boundary[i]][1] = (_enableTimeStep ? (c + 2) * _mag / _nt : _mag) * sin(mesh_points[bottom_boundary[i]][0]);
        }
    }
}

void springSmooth::outNewMesh(const double &current_t)
{
    unsigned int _points_num = grid.returnPointsNum();
    unsigned int _elem_num = grid.returnElemNum();

    std::vector<std::vector<double>> mesh_points = grid.retrunMesh_points();
    std::vector<std::vector<unsigned int>> _ele_topo = grid.return_elem_topo();

    std::ofstream fout("tec.plt", std::ios_base::app);
    fout << " VARIABLES = \" X \", \" Y \" " << std::endl;
    fout << "ZONE T=\"dynamicMesh\", NODES=" << _points_num << ", ELEMENTS=" << _elem_num << ", DATAPACKING=POINT, ZONETYPE=FEQuadrilateral,"
         << "SOLUTIONTIME=" << current_t << std::endl;

    for (int i = 0; i < _points_num; i++)
    {
        for (int j = 0; j < 2; j++)
        {
            mesh_points[i][j] += d_current[i][j];
        }

        fout << std::setw(15) << mesh_points[i][0] << std::setw(15) << mesh_points[i][1] << std::endl;
    }

    for (int i = 0; i < _elem_num; i++)
    {
        fout << std::setw(7) << _ele_topo[i][0] + 1 << std::setw(7) << _ele_topo[i][1] + 1 << std::setw(7)
             << _ele_topo[i][2] + 1 << std::setw(7) << _ele_topo[i][3] + 1 << std::endl;
    }

    fout.close();
}
